#	library(foreign);library(msm);library(MASS);library(VGAM);library(distr);library(truncnorm); library(Rcpp); library(RcppArmadillo); library(inline); library(RcppEigen)
	expit<-function(x) exp(x)/(1+exp(x))

#############################
#############################
#############################

criterion<-function(M, missing.mat=NULL, gibbs=100, burnin=100, max.optim=50, thin=1, save.curr="UDV_curr", save.each=FALSE, thin.save=25, maxdim=NULL,ratio.bin=1,EM=FALSE,nullmod=FALSE,prior=c(1,1.78)){

cat("Step 1 of 4: Formatting Data","\n")

	if(EM==TRUE) {max.optim<-1e10;burnin=1}
	votes.mat<-M
	if(length(missing.mat)==0) missing.mat<-votes.mat*0
	empir.cdf=NULL
#source('FunctionsInternal_Count_EM.R')
	test.gamma.pois<-test.gamma.pois_EM
	make.Z<-make.Z_EM
	update_UDV<-update_UDV_EM

#votes.mat=dtm2;missing.mat=dtm2*0;max.loop=500; burnin=50; thin=1;empir.cdf=NULL;cutoff.seq=seq.vec;save.curr="UDV_curr";save.each=FALSE;thin.save=25;max.optim=300

		dubcent<-function(X){
		colmean.vec<-colMeans(X,na.rm=T)
		rowmean.vec<-rowMeans(X,na.rm=T)
		
		for(i in 1:dim(X)[1]) for (j in 1:dim(X)[2]) 
		X[i,j]<-X[i,j]-rowmean.vec[i]-colmean.vec[j]
		
		X-mean(X)
		}

#if(cutoff.seq==0){
dtm2<-votes.mat
unique.dtm<-sort(unique(as.vector(dtm2)))

seq.vec<-c(-1:(max(dtm2)+2))
for(i in 2:length(seq.vec)) ifelse(seq.vec[i]%in%unique.dtm, seq.vec[i]<-seq.vec[i-1]+1,seq.vec[i]<-seq.vec[i-1] )

cutoff.seq<-seq.vec
rm(seq.vec)
#}




	n<-dim(votes.mat)[1]
	k<-dim(votes.mat)[2]
	##For only fitting first few dimensions; not in this version
	if(length(maxdim)==0) maxdim<-min(n,k) else maxdim<-min(n,k,maxdim)
	
	
	
	votes.mat<-as.matrix(votes.mat)
	
	row.type<-apply(votes.mat,1,FUN=function(x)  (max(x)<=1))
	row.type[row.type==TRUE]<-"bin"
	row.type[row.type==FALSE]<-"count"
	row.type.ord<-apply(votes.mat,1,FUN=function(x)  (min(x)<0))
	row.type[row.type.ord==TRUE]<-"ord"

tau.ord<-1
if(min(votes.mat)<0) {
	ord.all<-(votes.mat[row.type=="ord",])[missing.mat[row.type=="ord"]==0]
	tau.ord<-qnorm(mean(ord.all>-1))-qnorm(mean(ord.all==0))
}

	cat('Distribution of observed data\n')
	print(table(row.type))
	sigma<-1

##Initialize D, lambda, Theta, Z
	votes.mat.init<-t(apply(votes.mat,1,FUN=function(x) x/sd(x)))
	votes.mat.init[!is.finite(votes.mat.init)]<-0

	mu.r<-rowMeans(sigma*votes.mat.init)
	mu.c<-colMeans(sigma*votes.mat.init)
	mu.grand<- mean(sigma*votes.mat.init)
	
	ones.r<-rep(1,k)
	ones.c<-rep(1,n)
	
	
	svd.mat<- sigma*votes.mat.init - ones.c%*%t(mu.c)-mu.r%*%t(ones.r)+mu.grand
	svd.mat<- svd.mat-ones.c%*%t(colMeans(svd.mat))-rowMeans(svd.mat)%*%t(ones.r)+mean(svd.mat)
	
	
	Dtau<-rep(1, dim(votes.mat)[1])
	Theta.last<-Z.next<-Z.last<-0*votes.mat
	Theta.last<-Z.last<-U.last<-0*votes.mat

##Containers
	U.run<-U.run.2<-matrix(NA, nrow=dim(votes.mat)[1], ncol=gibbs+burnin)
	V.run<-V.run.2<-matrix(NA, nrow=dim(votes.mat)[2], ncol=gibbs+burnin)
	D.run<-D.mean.run<-matrix(NA, nrow=maxdim, ncol=gibbs+burnin)
	accept.run<-NULL
	lambda.run<-NULL
	waicsparse.run<-waicmat.run<-NULL
	dev.run<-matrix(NA, nrow=2, ncol=gibbs+burnin)
	V.average<-U.average<-0
	Z.run<-Z.mode.run<-0
	if(sum(row.type=="bin")>2){
	Theta.last[row.type=="bin",][votes.mat[row.type=="bin",]==1]<-dnorm(0)/pnorm(.5)
	Theta.last[row.type=="bin",][votes.mat[row.type=="bin",]==0]<- -dnorm(0)/pnorm(.5)
	Theta.last[row.type=="bin",][missing.mat[row.type=="bin",]==1]<- 0}
		dubcent<-function(X){
		colmean.vec<-colMeans(X,na.rm=T)
		rowmean.vec<-rowMeans(X,na.rm=T)
		
		for(i in 1:dim(X)[1]) for (j in 1:dim(X)[2]) 
		X[i,j]<-X[i,j]-rowmean.vec[i]-colmean.vec[j]
		
		X-mean(X)
		}

	if(sum(row.type=="count")>2){
	Theta.last[row.type=="count",]<- (log(1*(votes.mat[row.type=="count",]==0)+votes.mat[row.type=="count",]))
	#Theta.last[row.type=="count",]<-(Theta.last[row.type=="count",]-mean(Theta.last[row.type=="count",]))/sd(Theta.last[row.type=="count",])
	Theta.last[row.type=="count",][missing.mat[row.type=="count",]==1]<- 0
	#Theta.last[row.type=="count",]<- 0
	}


	if(sum(row.type=="ord")>2){
	Theta.last[row.type=="ord",][votes.mat[row.type=="ord",]==-2]<-dnorm(0)/pnorm(.5)
	Theta.last[row.type=="ord",][votes.mat[row.type=="ord",]==0]<- -dnorm(0)/pnorm(.5)
	Theta.last[row.type=="ord",][votes.mat[row.type=="ord",]==-1]<- 0

	Theta.last[row.type=="ord",][missing.mat[row.type=="ord",]==1]<- 0}


	
	svd0<-svd(svd.mat,nu=maxdim,nv=maxdim)
	UDV.all<-NULL
	UDV.all$V.next<-svd0$v
	lambda.shrink<-lambda.lasso<-rep(median((2/(svd0$d))^.5),maxdim)[1]


cat("Step 2 of 4: Beginning Burnin Period","\n")


	params.init<-c(0,0)#log(c(mean(votes.mat)))
	proposal.sd<-0.1
	step.size<-scale.sd<-1
	waic.count<-1
	taus.run<-0
#############################
#############################
#############################		
###########Start loop here

for(i.gibbs in 1:(gibbs+burnin)){
	#Switch off EM
	if(i.gibbs==max.optim) {
		test.gamma.pois<-test.gamma.pois_gibbs
		make.Z<-make.Z_gibbs
		update_UDV<-update_UDV_gibbs	
			
#		source('FunctionsInternal_Count.R')
	}

	#Step 1: Make Z-star
	Z.call<-make.Z(Theta.last.0=Theta.last, votes.mat.0=votes.mat, row.type.0=row.type,  n0=n ,k0=k, iter.curr=i.gibbs,empir=empir.cdf,cutoff.seq.0=cutoff.seq,missing.mat.0=missing.mat,lambda.lasso=lambda.lasso,params=params.init,proposal.sd=proposal.sd,scale.sd=scale.sd,max.optim=max.optim,step.size=step.size,maxdim.0=maxdim,tau.ord.0=tau.ord,test.gamma.pois.0=test.gamma.pois)


	##Capture WAIC
	if(i.gibbs==max(max.optim,burnin)) {
		waic.mat<-Z.call$waicmat
		taus.run<-taus.run+Z.call$tau.ord
		}
	if(i.gibbs>max(max.optim,burnin)){
		waic.mat[,2:4]<-waic.mat[,2:4]+Z.call$waicmat[,2:4]
		taus.run<-taus.run+Z.call$tau.ord
		waic.count<-waic.count+1
	}
	
	if(EM==TRUE) {
		taus.mean<-taus.run<-Z.call$tau.ord
		}
		
	Z.next<-Z.call$Z
	params.init<-Z.call$params.init
	accept.run[i.gibbs]<-Z.call$accept

	#dev.run[i.gibbs]<-Z.call$prob
	dev.run[,i.gibbs]<-Z.call$dev
	converge<-FALSE
	if(EM==TRUE&i.gibbs>20){
		if(max(abs(dev.run[(i.gibbs-19):i.gibbs]-dev.run[i.gibbs]))<0.1  ) {
			print("EM Converged")
			converge<-TRUE
			}
		
	}
	if(converge) break
	params.init<-Z.call$params
	proposal.sd<-Z.call$proposal.sd
	step.size<-Z.call$step.size
	tau.ord<-Z.call$tau.ord
	#Step 2: Update UDV and components therein
	UDV.all<-	update_UDV(
	Z.next.0=Z.next,
	k0=k, n0=n, lambda.lasso.0=lambda.lasso,lambda.shrink.0=lambda.shrink,
	votes.mat.0=votes.mat, Dtau.0=Dtau, iter.curr=i.gibbs, row.type.0=row.type, missing.mat.0=missing.mat,maxdim.0=maxdim, V.last = UDV.all$V.next,ratio.bin.0=ratio.bin,prior.0=prior
	)
	
	last.Z<-Z.call
	last.UDV<-UDV.all
	#print(tau.ord)
	
	
	#Change--update stepsize during burnin only
	if(i.gibbs>(max.optim+50)&i.gibbs%%50==1){
	#if(i.gibbs>(max.optim+50)&i.gibbs%%50==1 & i.gibbs<burnin ){

		seq.check<-(i.gibbs-1):(i.gibbs-20)
		#print(accept.run)
		if(mean(accept.run[(i.gibbs-1):(i.gibbs-20)])>.9 ) step.size<-step.size*1.25
		if(mean(accept.run[(i.gibbs-1):(i.gibbs-20)])<.1 ) step.size<-step.size*.8
		#print("acceptance percentage")
		#print(mean(accept.run[(i.gibbs-1):(i.gibbs-20)]))
	}
	

	if(save.each){
          file.temp<-paste(save.curr,i.gibbs)
          if(i.gibbs%%thin.save ==1) {save(UDV.all,file=file.temp)}
	} else {
          ## file.temp<-paste("/scratch/insong/sfa_output/", save.curr,i.gibbs, sep="")
          file.temp <- paste(save.curr)
          save(UDV.all,file=file.temp)
	}
	#Update inputs to UDV.all
	Dtau=UDV.all$Dtau
	if(nullmod) Dtau<-Dtau*0+1e-10
	lambda.lasso<-UDV.all$lambda.lasso
	lambda.shrink<-UDV.all$lambda.shrink
	Theta.last<-UDV.all$Theta.last

	
	
	if(i.gibbs==1) cat("    0% of Burnin Completed:", i.gibbs,"/",burnin,"\n")	
	if(i.gibbs==floor(burnin/4)) cat("   25% of Burnin Completed:", i.gibbs,"/",burnin,"\n")
	if(i.gibbs==floor(burnin/2)) cat("   50% of Burnin Completed:", i.gibbs,"/",burnin,"\n")
	if(i.gibbs==floor(3*burnin/4)) cat("   75% of Burnin Completed:", i.gibbs,"/",burnin,"\n")
	if(i.gibbs==floor(burnin)) cat("  100% of Burnin Completed:", i.gibbs,"/",burnin,"\n")
	if(i.gibbs == burnin+1) cat("Step 3 of 4: Burnin Completed; Saving Samples \n")
	if(i.gibbs==floor(burnin+1)) {
		cat("    0% of Samples Gathered:", i.gibbs-burnin,"/",gibbs,"\n")
		cat("Current dimensionality estimate:\n")
		print(UDV.all$D.trunc[1:maxdim])
		}
	if(i.gibbs==floor(burnin+(gibbs)/4)) {
		cat("   25% of Samples Gathered:", i.gibbs-burnin,"/",gibbs,"\n")
		cat("Current dimensionality estimate:\n")
		print(UDV.all$D.trunc[1:maxdim])
		}
	if(i.gibbs==floor(burnin+(gibbs)/2)) {
		cat("   50% of Samples Gathered:", i.gibbs-burnin,"/",gibbs,"\n")
		cat("Current dimensionality estimate:\n")
		print(UDV.all$D.trunc[1:maxdim])
		}
	if(i.gibbs==floor(burnin+3*(gibbs)/4)) {
		cat("   75% of Samples Gathered:", i.gibbs-burnin,"/",gibbs,"\n")
		cat("Current dimensionality estimate:\n")
		print(UDV.all$D.trunc[1:maxdim])
		}
	if(i.gibbs==floor(burnin+gibbs)) {
		cat("  100% of Samples Gathered:", i.gibbs-burnin,"/",burnin,"\n")
		cat("Current dimensionality estimate:\n")
		print(UDV.all$D.trunc[1:maxdim])
		}

##Gather results
	lambda.run[i.gibbs]<-UDV.all$lambda.lasso[1]
	U.run[,i.gibbs]<-UDV.all$U.next[,1]
	U.run.2[,i.gibbs]<-UDV.all$U.next[,2]
	V.run[,i.gibbs]<-UDV.all$V.next[,1]
	V.run.2[,i.gibbs]<-UDV.all$V.next[,2]
	D.run[,i.gibbs]<-UDV.all$D.trunc[1:maxdim]
	D.mean.run[,i.gibbs]<-UDV.all$D.post[1:maxdim]
	

	
	post.process.inner<-function(x,y){
			for(i.cor in 1:ncol(x)) y[,i.cor]<-y[,i.cor]*sign(cov(x[,i.cor],y[,i.cor]))
			return(y)
	}
	if(i.gibbs==burnin){
		U.order<-UDV.all$U.next
		V.order<-UDV.all$V.next
		U.average<-V.average<-0
	}

	if(i.gibbs>burnin & i.gibbs %% thin == 0) {
		Z.run<- Z.run + Theta.last
		Z.mode.run<- Z.mode.run + UDV.all$Theta.mode
	
		##WAIC off mode
		Z.one<- UDV.all$Theta.mode[row.type=="bin",][votes.mat[row.type=="bin",]==1]
		Z.zero<- UDV.all$Theta.mode[row.type=="bin",][votes.mat[row.type=="bin",]==0]
		Z.one[abs(Z.one)>38]<-38*sign(Z.one[abs(Z.one)>38])
		Z.zero[abs(Z.zero)>38]<-38*sign(Z.zero[abs(Z.zero)>38])
		
		prob.vec<-as.vector(c(pnorm(-Z.zero,log=TRUE ),pnorm(Z.one,log=TRUE )))
		voteprob.mat<-data.frame("bin",exp(prob.vec),prob.vec,prob.vec^2)
		
		names(voteprob.mat)<-c("type","p","lp","psq")
		
		vote.sparsedev<--2*sum(prob.vec,na.rm=T)

		pois.obj<-test.gamma.pois(params.init,votes.mat.0=votes.mat[row.type=="count",], Theta.last.0= UDV.all$Theta.mode[row.type=="count",],cutoff.seq=cutoff.seq)
		taus<-pois.obj$tau
		#print(pois.obj$data)
		word.sparsedev<-pois.obj$data
		prob.vec<-as.vector(pois.obj$probvec)
		
		wordprob.mat<-data.frame("count",exp(prob.vec),prob.vec,prob.vec^2)
		names(wordprob.mat)<-c("type","p","lp","psq")
		
		waic.temp<-rbind(voteprob.mat,wordprob.mat)
		if(i.gibbs==burnin+1){
			waicsparse.run<-waic.temp
			waicsparse.run[,2:4]<-0
		}
		
		waicsparse.run[,1]<-waic.temp[,1]
		waicsparse.run[,2:4]<-waicsparse.run[,2:4]+waic.temp[,2:4]/gibbs


		
		U.average<-U.average+post.process.inner(U.order,UDV.all$U.next)
		V.average<-V.average+post.process.inner(V.order,UDV.all$V.next)
		}
	#if(i.gibbs %%100 ==0) print(i.gibbs)
}#Closes out for loop
	cat("Step 4 of 4: Gathering Output\n")
	
	
	##Gather waic
	if(EM==FALSE)		{
		waic.mat[,2:4]<-waic.mat[,2:4]/waic.count
		taus.mean<-taus.run/waic.count
		WAIC <- -2*sum(log(waic.mat[,2]),na.rm=TRUE)+2*(sum(waic.mat[,4]-waic.mat[,3]^2))
		taus.mean<-taus.run/waic.count
		
		sparseWAIC<--2*sum(log(waicsparse.run[,2]),na.rm=TRUE)+2*(sum(waicsparse.run[,4]-waicsparse.run[,3]^2))
		}else{
		waic.mat<-NULL
		WAIC<-NULL
			}

				
	post.process<-function(x){apply(x,2,FUN=function(x2) x2*sign(cor(x2,x[,1])) ) }
	V.run<-post.process(V.run[,-c(1:burnin)])
	V.run.2<-post.process(V.run.2[,-c(1:burnin)])
	U.run <-post.process(U.run[,-c(1:burnin)])
	U.run.2 <-post.process(U.run.2[,-c(1:burnin)])



##Summary.stats

out<-list("dim.sparse"=D.run[,-c(1:burnin)], "dim.mean"=D.mean.run[,-c(1:burnin)], "rowdim1"=U.run, "rowdim2"=U.run.2, 
"coldim1"=V.run, "coldim2"=V.run.2, 
"lambda.lasso"=lambda.run[-c(1:burnin)], "Z" = Z.run/gibbs, "Z.mode"= Z.mode.run/gibbs,
"rowdims.all"=U.average/gibbs, "coldims.all"=V.average/gibbs,"dev"=dev.run,"waic.out"=waic.mat,"taus"=taus.mean,"WAIC"=WAIC,"sparseWAIC"=sparseWAIC,"waic.sparsemat"=waicsparse.run)
class(out)<-"sfa"
invisible(out)
}


################################
################################
##Three functions in here:
#### 1) test.gamma for finding the ML estimate for the grand lambda
#### 2) make.Z for converting a matrix to the Z-scale
#### 3) update.UDV for updating the ideal point estimates
#### *) Small additonal ones at the end 
################################
################################

################################
################################
##1) Finding gamma
#test.gamma.pois_EM<-function(gamma.try,Theta.last.0=Theta.last[row.type=="count",],votes.mat.0=votes.mat[row.type=="count",],emp.cdf.0,cutoff.seq=NULL){
test.gamma.pois_EM<-function(gamma.try,Theta.last.0,votes.mat.0,emp.cdf.0,cutoff.seq=NULL){	
	gamma.try<-exp(gamma.try)
	votes.mat<-votes.mat.0
	taus.try<-NULL
	count.seq<-cutoff.seq
	if(length(cutoff.seq)==0) count.seq<-seq(-1,max(votes.mat)+2,1)
	taus.try<-count.seq*0

	analytic.cdf<-count.seq
	a.int<-qnorm(mean(votes.mat==0))
	taus.try[count.seq>=0]<-(a.int+gamma.try[1]*count.seq[count.seq>=0]^(gamma.try[2]))
	taus.try[count.seq<0]<- -Inf
	taus.try<-sort(taus.try)
	taus.try[1]<--Inf	
	find.pnorm<-function(x){
		a<-x[1]
		b<-x[2]
		coefs<-c( -1.82517672, 0.51283415, -0.81377290, -0.02699400, -0.49642787, -0.33379312, -0.24176661, 0.03776971)
		x<-c(1, a, b, a^2,b^2, log(abs(a-b)),log(abs(a-b))^2,a*b)
		(sum(x*coefs))
	}
	lik<-((pnorm(taus.try[votes.mat+2 ]-Theta.last.0)-pnorm(taus.try[votes.mat+1 ]-Theta.last.0)) )
	log.lik<-log(lik)
	
	which.zero<-which(lik==0)
	a0<-taus.try[votes.mat[which.zero]+2]-Theta.last.0[which.zero]
	b0<-taus.try[votes.mat[which.zero]+1]-Theta.last.0[which.zero]
	log.lik[which.zero]<-apply(cbind(a0,b0),1,find.pnorm)
	log.lik[is.infinite(lik)&lik>0]<-max(log.lik[is.finite(lik)])
	log.lik[is.infinite(lik)&lik<0]<-min(log.lik[is.finite(lik)])

	thresh<-min(-1e30, min(log.lik[is.finite(log.lik)],na.rm=TRUE))
	log.lik[log.lik<thresh]<-thresh
	data.dev<--2*sum(log.lik,na.rm=TRUE)
	dev.out<-data.dev+sum(log(gamma.try)^2)
	
	
	return(list("deviance"=dev.out,"tau"=taus.try,"data.dev"=data.dev))

}	



################################
################################
##2) Converting a matrix to the Z scale
	make.Z_EM<-function(
			Theta.last.0=Theta.last, 
			votes.mat.0=votes.mat,row.type.0=row.type, 
			n0=n, k0=k, params=NULL, iter.curr=0,empir=NULL,cutoff.seq.0=NULL,missing.mat.0=NULL,lambda.lasso,proposal.sd,scale.sd, max.optim,step.size,maxdim.0,tau.ord.0,test.gamma.pois.0
			){
		
		tau.ord<-tau.ord.0
		Theta.last<-Theta.last.0;votes.mat<-votes.mat.0;
		row.type<-row.type.0; n<-n0; k<-k0
		cutoff.seq<-cutoff.seq.0
		sigma<-1
		row.type<-row.type.0; votes.mat<-votes.mat.0
		Z.next<-matrix(NA,nrow=n,ncol=k)
		missing.mat<-missing.mat.0
		i.gibbs<-iter.curr
		maxdim<-maxdim.0
		#print(i.gibbs)
		#print(iter.curr)
	test.gamma.pois<-test.gamma.pois.0
	
	
		##Initialize deviance
		word.dev<-vote.dev<-0


	estep.bin<-function(means,a,b){
		means+(dnorm(a-means)-dnorm(b-means))/(pnorm(b-means)-pnorm(a-	means))
		}
		
		
	
			if(sum(row.type=="bin")>0){
		toplimit.mat<-bottomlimit.mat<-votes.mat[row.type=="bin",]*0
		toplimit.mat[votes.mat[row.type=="bin",]==1]<-Inf
		toplimit.mat[votes.mat[row.type=="bin",]==0]<-0
		toplimit.mat[votes.mat[row.type=="bin",]==0.5]<-Inf
		bottomlimit.mat[votes.mat[row.type=="bin",]==1]<-0
		bottomlimit.mat[votes.mat[row.type=="bin",]==0]<--Inf
		bottomlimit.mat[votes.mat[row.type=="bin",]==0.5]<- -Inf


		Z.next[row.type=="bin",]<-
estep.bin(means=Theta.last[row.type=="bin",],a=bottomlimit.mat,b=toplimit.mat)
		Z.next[row.type=="bin",][is.na(Z.next[row.type=="bin",])]<-0
		Z.next[row.type=="bin",][is.infinite(Z.next[row.type=="bin",])]<-0

		
		pars.max<-1
		accept.out<-prob.accept<-1
		Z.one<-Z.next[row.type=="bin",][votes.mat[row.type=="bin",]==1]
		Z.zero<-Z.next[row.type=="bin",][votes.mat[row.type=="bin",]==0]
		Z.one[abs(Z.one)>38]<-38*sign(Z.one[abs(Z.one)>38])
		Z.zero[abs(Z.zero)>38]<-38*sign(Z.zero[abs(Z.zero)>38])
		
		word.dev<--2*(sum(pnorm(-Z.zero,log=TRUE ),na.rm=TRUE)+sum(pnorm(Z.one,log=TRUE ),na.rm=TRUE))
		
		
		
		}
		
		
			if(sum(row.type=="ord")>0){
				
			#Z.next[row.type=="ord",][votes.mat[row.type=="ord",]==-2]<-rtnorm(sum(votes.mat[row.type=="ord",]==-2), mean=sigma^.5*Theta.last[row.type=="ord",][votes.mat[row.type=="ord",]==-2], lower=tau.ord, sd=sigma^.5)
			#Z.next[row.type=="ord",][votes.mat[row.type=="ord",]==-1]<-rtnorm(sum(votes.mat[row.type=="ord",]==-1), mean=sigma^.5*Theta.last[row.type=="ord",][votes.mat[row.type=="ord",]==-1], lower=0,upper=tau.ord, sd=sigma^.5)
		#Z.next[row.type=="ord",][votes.mat[row.type=="ord",]==0]<-rtnorm(sum(votes.mat[row.type=="ord",]==0), mean=sigma^.5*Theta.last[row.type=="ord",][votes.mat[row.type=="ord",]==0], upper=0 , sd=sigma^.5)
		#Z.next[row.type=="ord",][missing.mat[row.type=="ord",]==1]<-rnorm(sum(missing.mat[row.type=="ord",]==1), mean=sigma^.5*Theta.last[row.type=="ord",][missing.mat[row.type=="ord",]==1], sd=sigma^.5)
		
		
		lower.mat<-upper.mat<-matrix(0,nrow=nrow(votes.mat[row.type=="ord",]), ncol=ncol(votes.mat))
		##Top category
		lower.mat[votes.mat[row.type=="ord",]==-2]<- tau.ord
		upper.mat[votes.mat[row.type=="ord",]==-2]<- Inf
		##Middle category
		lower.mat[votes.mat[row.type=="ord",]==-1]<- 0
		upper.mat[votes.mat[row.type=="ord",]==-1]<- tau.ord
		##Lower category
		lower.mat[votes.mat[row.type=="ord",]==0]<- -Inf
		upper.mat[votes.mat[row.type=="ord",]==0]<- 0
		#Missing data
		lower.mat[missing.mat[row.type=="ord",]==1]<- -Inf
		upper.mat[missing.mat[row.type=="ord",]==1]<- Inf

	
		Z.next.temp<-estep.bin(means=Theta.last[row.type=="ord",],a=lower.mat,b=upper.mat)
		which.change<-!is.finite(Z.next.temp)	
		Z.next.temp[which.change]<-Theta.last[row.type=="ord",][which.change]
		Z.next[row.type=="ord",]<-Z.next.temp



		tau.min<-max(Z.next[row.type=="ord",][votes.mat[row.type=="ord",]==-1&missing.mat[row.type=="ord"]==0])
		tau.max<-min(Z.next[row.type=="ord",][votes.mat[row.type=="ord",]==-2&missing.mat[row.type=="ord"]==0])
		tau.min<-max(tau.min,0.001)		
		tau.samp<-sort(c(tau.ord,runif(100,min(tau.min,tau.max),max(tau.min,tau.max))))
		#which.tau<-which(tau.samp==tau.ord)
		#tau.ord<-tau.samp[102-which.tau]
		#print(c(tau.min,tau.max))
		tau.ord<-runif(1,min(tau.min,tau.max),max(tau.min,tau.max))
		pars.max<-1
		accept.out<-prob.accept<-1

				}

		if(sum(row.type=="count")>0){
		#begin type = "pois"


	num.sample.count<-sum(row.type=="count")*k
	accept.out<-prob.accept<-NA

	dev.est<-function(x) test.gamma.pois_EM(x,votes.mat.0=votes.mat[row.type=="count",], Theta.last.0=Theta.last[row.type=="count",],cutoff.seq=cutoff.seq)$de
	dev.est1<-function(x) test.gamma.pois_EM(c(x,params[2]),votes.mat.0=votes.mat[row.type=="count",], Theta.last.0=Theta.last[row.type=="count",],cutoff.seq=cutoff.seq)$de
	dev.est2<-function(x) test.gamma.pois_EM(c(params[1],x),votes.mat.0=votes.mat[row.type=="count",], Theta.last.0=Theta.last[row.type=="count",],cutoff.seq=cutoff.seq)$de
	
		params<-round(params,dig=2)

	#range.opt<-c(params-1,params+1)
	range.opt<-rbind(params-1,params+1)
	if(iter.curr>10) range.opt<-rbind(params-.25,params+.25)

	##Change -- code block below did maximization. Let's do E step
	if(TRUE){
	if(iter.curr%%1==0|iter.curr<3){
	gamma.opt.1<-optimize(dev.est1,
			lower=range.opt[1,1],upper=range.opt[2,1],tol=0.01)
	params[1]<-gamma.opt.1$minimum
	gamma.opt.2<-optimize(dev.est2,
		lower=range.opt[1],upper=range.opt[2],tol=0.01)
	params[2]<-gamma.opt.2$minimum
	gamma.opt<-list(gamma.opt.1,gamma.opt.2)
	"M Step"
	#print(gamma.opt.2)

	}
	}
	
	
	params<-round(params,dig=2)
	
	##Change -- new E step is here below
	if(FALSE){
	params<-NULL
	range.opt<-round(range.opt,dig=3)
	
	##First param
	pars1<-seq(range.opt[1,1],range.opt[1,2],length=20)
	dev1<-sapply(pars1,dev.est1)
	probs1<-dev1/-2
	probs1<-exp(probs1-max(probs1))
	
	params[1]<-sum(pars1*probs1)/sum(probs1)

	##First param
	pars2<-seq(range.opt[2,1],range.opt[2,2],length=20)
	dev2<-sapply(pars2,dev.est2)
	probs2<-dev2/-2
	probs2<-exp(probs2-max(probs2))
	
	params[2]<-sum(pars2*probs2)/sum(probs2)
	}
	
	params<-round(params,dig=2)
	
	
	gamma.next<-pars.max<-params

		#taus<-test.gamma.pois_EM(gamma.next,votes.mat.0=votes.mat[row.type=="count",], Theta.last.0=Theta.last[row.type=="count",],cutoff.seq=cutoff.seq)$tau

		pois.obj<-test.gamma.pois(gamma.next,votes.mat.0=votes.mat[row.type=="count",], Theta.last.0=Theta.last[row.type=="count",],cutoff.seq=cutoff.seq)
		taus<-pois.obj$tau
		#print()
		word.dev<-pois.obj$data
		taus<-sort(taus)
		taus[1]<--Inf
		#print("Range of gamma")
		#print(range.opt)
		#print(gamma.next)
		#print(range(taus))

#function(means,a,b){
#	means+(dnorm(a-means)-dnorm(b-means))/(pnorm(b-means)-pnorm(a-means))
	#}		
		
		lower.mat<-matrix(taus[votes.mat[row.type=="count",]+1 ],nrow=nrow(votes.mat[row.type=="count",]))
		upper.mat<-matrix(taus[votes.mat[row.type=="count",]+2 ],nrow=nrow(votes.mat[row.type=="count",]))

	

		Z.next.temp<-estep.bin(means=Theta.last[row.type=="count"],a=lower.mat,b=upper.mat)

		which.change<-is.na(Z.next.temp)|is.infinite(Z.next.temp)

		Z.next.temp[which.change]<-
		(Theta.last[row.type=="count",][which.change] > upper.mat[which.change])*
		upper.mat[which.change]+
		(Theta.last[row.type=="count",][which.change] < lower.mat[which.change])*lower.mat[which.change]
		
		which.change<-is.na(Z.next.temp)|is.infinite(Z.next.temp)
		Z.next.temp[which.change]<-Theta.last[which.change]
		
				
		Z.next[row.type=="count",]<-Z.next.temp
		}
		
		dev.all<-word.dev+vote.dev
		
	return(list("Z.next"=Z.next,"params"=(pars.max),"accept"=accept.out,"prob"=prob.accept,"proposal.sd"=proposal.sd,"step.size"=step.size,"tau.ord"=tau.ord,"dev"=c(word.dev,vote.dev)))
	}##Closes out make.Z function


################################
################################
##2) Converting a matrix to the Z scale
	update_UDV_EM<-function(
	Z.next.0=Z.next,
	k0=k, n0=n, lambda.lasso.0=lambda.lasso,lambda.shrink.0=lambda.shrink,
	Dtau.0=Dtau,
	votes.mat.0=votes.mat, iter.curr=0,row.type.0,missing.mat.0=missing.mat,maxdim.0,
	V.last,ratio.bin.0, prior.0
	){

	missing.mat<-missing.mat.0
	Dtau<-Dtau.0;
	Z.next<-Z.next.0;
	votes.mat<-votes.mat.0;
	n<-n0; k<-k0; 
	lambda.lasso<-lambda.lasso.0
	row.type <- row.type.0
	maxdim<-maxdim.0
	prior<-prior.0

	#Declare some vectors
	sigma<-1
	ones.r<-rep(1,k0)
	ones.c<-rep(1,n0)

	#Update intercepts
	mu.r<-rowMeans(Z.next)#-ones.c%*%t(mu.c)-Theta.last.0+mu.grand)
	mu.r<-mu.r*n/(n+1)#+rnorm(length(mu.r),sd=1/k)

	mu.c<-colMeans(Z.next)#-mu.r%*%t(ones.r)-Theta.last.0+mu.grand)
	#mu.c<-mu.c*k/(k+1)#+rnorm(length(mu.c),sd=1/n)
	mu.grand<-mean(Z.next)
	mu.grand<-mu.grand*(n*k)/(n*k+1)#+rnorm(1,sd=1/(n*k))
	
	mean.mat<-ones.c%*%t(mu.c)+mu.r%*%t(ones.r)-mu.grand
	
	if(length(unique(row.type))>1){
		is.bin<-sum(row.type=="bin")>0
		is.count<-sum(row.type=="count")>0
		is.ord<-sum(row.type=="ord")>0
		
		#print("Two Means Being Used")
		mean.c.mat<-matrix(NA,nrow=n,ncol=k)
		if(is.bin) mu.c.bin<-colMeans(Z.next[row.type=="bin",])
		if(is.count) mu.c.count<-colMeans(Z.next[row.type=="count",])
		if(is.ord) mu.c.ord<-colMeans(Z.next[row.type=="ord",])

		if(is.bin) mu.c.bin<-(mu.c.bin*k)/(k+1)#+rnorm(length(mu.c.bin),sd=1/sum(row.type=="bin"))
		if(is.count) mu.c.count<-mu.c.count*k/(k+1)#+rnorm(length(mu.c.count),sd=1/sum(row.type=="count"))
		if(is.ord) mu.c.ord<-mu.c.ord*k/(k+1)#+rnorm(length(mu.c.count),sd=1/sum(row.type=="count"))

		if(is.bin) mean.c.mat[row.type=="bin",]<-ones.c[row.type=="bin"]%*%t(mu.c.bin)
		if(is.count) mean.c.mat[row.type=="count",]<-ones.c[row.type=="count"]%*%t(mu.c.count)
		if(is.ord) mean.c.mat[row.type=="ord",]<-ones.c[row.type=="ord"]%*%t(mu.c.ord)

	
		mean.grand.mat<-matrix(NA,nrow=n,ncol=k)
		if(is.bin) mean.grand.mat[row.type=="bin",]<-mean(Z.next[row.type=="bin",])* (sum(row.type=="bin")*k)/(sum(row.type=="bin")*k+1)#+rnorm(1,sd=1/(sum(row.type=="bin")*k))
		if(is.count) mean.grand.mat[row.type=="count",]<-mean(Z.next[row.type=="count",])* (sum(row.type=="count")*k)/(sum(row.type=="count")*k+1)#+rnorm(1,sd=1/(sum(row.type=="bin")*k))#+rnorm(1,sd=1/(sum(row.type=="count")*k))
		if(is.ord) mean.grand.mat[row.type=="ord",]<-mean(Z.next[row.type=="ord",])* (sum(row.type=="count")*k)/(sum(row.type=="ord")*k+1)#+rnorm(1,sd=1/(sum(row.type=="bin")*k))#+rnorm(1,sd=1/(sum(row.type=="count")*k))


	mean.mat<-mean.c.mat+mu.r%*%t(ones.r)-mean.grand.mat
	
	}

	Z.starstar<-svd.mat<- Z.next- mean.mat




	

	#Take svd, give each column an sd of 1 (rather than norm of 1)
	num.zeroes<-1#colMeans(votes.mat!=0)
	svd.mat.0<-svd.mat

	
	#save(svd.mat,file="svd.mat")
	
	svd.mat[is.na(svd.mat)]<-0
	
	wts.dum<-rep(1,nrow(svd.mat))
	wts.dum[row.type=="bin"]<-1/sum(1-missing.mat[row.type=="bin",])^.5*ratio.bin.0
	wts.dum[row.type=="count"]<-1/sum(1-missing.mat[row.type=="count",])^.5
	wts.dum[row.type=="ord"]<-1/sum(1-missing.mat[row.type=="ord",])^.5


	wts.dum<-wts.dum/mean(wts.dum)

	#svd.dum<-irlba(svd.mat*wts.dum,nu=maxdim,nv=maxdim,V=V.last)
	##Change
	svd.dum<-svd(svd.mat*wts.dum,nu=maxdim,nv=maxdim)
	#svd.dum<-rsvd(svd.mat*wts.dum,k=maxdim)

	svd.dum$u[is.na(svd.dum$u)|is.infinite(svd.dum$u)]<-0
	svd.dum$d[is.na(svd.dum$d)|is.infinite(svd.dum$d)]<-0
	svd.dum$v[is.na(svd.dum$v)|is.infinite(svd.dum$v)]<-0


	svd.dum$d<-(t(svd.dum$u)%*%svd.mat%*%svd.dum$v)
	which.rows<-which(rowMeans(svd.dum$d^2)^.5<1e-4)
	which.cols<-which(colMeans(svd.dum$d^2)^.5<1e-4)
	##Change--put zeroes in for missing
	svd.dum$d[which.rows,]<-rnorm(length(svd.dum$d[which.rows,]),sd=.001)*0
	svd.dum$d[which.cols,]<-rnorm(length(svd.dum$d[which.cols,]),sd=.001)*0
	#Change
	svd2<-svd(svd.dum$d,nu=maxdim,nv=maxdim)	
	#svd2<-rsvd(svd.dum$d,k=maxdim)
	#print(svd.dum$d)
	#print(maxdim)
	#print(maxdim-2)
	#svd2<-irlba(svd.dum$d,nu=maxdim-5,nv=maxdim-5)
	svd2$u[is.na(svd2$u)|is.infinite(svd2$u)]<-0
	svd2$d[is.na(svd2$d)|is.infinite(svd2$d)]<-0
	svd2$v[is.na(svd2$v)|is.infinite(svd2$v)]<-0
	
	svd.dum$u<-svd.dum$u%*%(svd2$u)
	svd.dum$v<-svd.dum$v%*%(svd2$v)

	svd.dum$d<-svd2$d
	
	svd0<-svd.dum
	
	svd0$v<-t(t(svd0$u)%*%svd.mat.0)
	svd0$v<-apply(svd0$v,2,FUN=function(x) x/sum(x^2)^.5)
	svd0$d<-diag(t(svd0$u)%*%svd.mat.0%*%svd0$v)
	sort.ord<-sort(svd0$d,ind=T,decreasing=T)$ix
	svd0$u<-svd0$u[,sort.ord]
	svd0$v<-svd0$v[,sort.ord]
	svd0$d<-svd0$d[sort.ord]
	svd0$u<-svd0$u*(n-1)^.5
	svd0$v<-svd0$v*(k-1)^.5
	svd0$d<-svd0$d*((n-1)*(k-1))^-.5
	Theta.last.0<-svd.mat
	Theta.last<-Theta.last.0+(ones.c%*%t(mu.c)+mu.r%*%t(ones.r)-mu.grand)

	#Update d; follows from Blasso and DvD
	Y.tilde<-as.vector(svd.mat)
	if(n>length(Dtau)) Dtau[(length(Dtau)+1):n]<-1
	A<- (n*k)*diag(n)+diag(as.vector(Dtau^(-1)))
	gA<-A*0+NA
	gA[1:maxdim,1:maxdim]<-ginv(A[1:maxdim,1:maxdim])
	gA[is.na(gA)]<-0
	XprimeY<-sapply(1:maxdim, FUN=function(i, svd2=svd0, Z.use=svd.mat) sum(Z.use*(svd2$u[,i]%*%t(svd2$v[,i]))))
	if(length(XprimeY)<dim(gA)[1]) XprimeY[(length(XprimeY)+1) :dim(gA)[1]]<-0
	D.post.mean<- as.vector(gA%*%XprimeY)
	D.post.var.2<-gA

##Sample D and reconstruct theta, putting intercepts back in
	D.post<-rep(NA,length(D.post.mean))
	D.post[1:maxdim]<-D.post.mean[1:maxdim] + diag(D.post.var.2[1:maxdim,1:maxdim]) ^.5 #as.vector(mvrnorm(1, mu=D.post.mean[1:maxdim], D.post.var.2[1:maxdim,1:maxdim] ) )/sigma^.5
	D.post[is.na(D.post)]<-0
	abs.D.post<-abs(D.post)
	#abs.D.post.loop<-abs(mvrnorm(50, mu=D.post.mean[1:maxdim], D.post.var.2[1:maxdim,1:maxdim] ) )
	#print(abs.D.post.loop)

##Calculate MAP and mean estimate
	D.trunc<-pmax(svd0$d-lambda.lasso.0,0)
	
	U.last<-svd0$u
	V.last<-svd0$v


	U.mean.next<-0*U.last
	V.mean.next<-0*V.last
##Construct U and V
#prior var of 2

	#D.adj<-abs(D.post[1:maxdim])/svd0$d
	#U.last<-t(t(U.last)*(D.post[1:maxdim]^2/(D.post[1:maxdim]^2+1/(4*k))))
	#V.last<-t(t(V.last)*D.post[1:maxdim]^2/(D.post[1:maxdim]^2+1/(4*n)))
	U.last[!is.finite(U.last)]<-0
	V.last[!is.finite(V.last)]<-0

	U.next<-U.last
	V.next<-V.last

	Theta.last.0<-U.next%*%diag(D.post[1:maxdim])%*%t(V.next)
	Theta.last<-Theta.last.0+mean.mat
	Theta.mode<- U.next%*%diag(D.trunc[1:maxdim])%*%t(V.next)+mean.mat
	##Update muprime, invTau2, lambda.lasso
	muprime<-(abs(lambda.lasso*sqrt(sigma)))/abs.D.post#*colMeans(1/abs.D.post.loop)
	invTau2<-muprime#sapply(1:maxdim, FUN=function(i) rinv.gaussian(1, muprime[i], (lambda.lasso^2 ) ) )
	#invTau2<-matrix(NA,nrow=500,ncol=maxdim)
	invTau2<-sapply(1:maxdim, FUN=function(i) rinv.gaussian(500, muprime[i], (lambda.lasso^2 ) ) )
	Dtau<-colMeans(1/abs(invTau2))
	#lambda.lasso<-((maxdim)/(sum(Dtau[1:maxdim])/2))^.5
	#lambda.lasso<-(2/mean(Dtau))^.5
	#ran.gamma<-rgamma(1000, shape=maxdim+1  , rate=sum(Dtau[1:maxdim])/2+1.78  )
	#d1<-density(ran.gamma^.5,cut=0)
	#lambda.shrink<-lambda.lasso<-d1$x[d1$y==max(d1$y)]
	lambda.shrink<-lambda.lasso<-((maxdim-1+prior[1])/(sum(Dtau[1:maxdim])/2+prior[2]))^.5
	#lambda.shrink<-lambda.lasso<-rgamma(1, shape=maxdim+1  , rate=sum(Dtau[1:maxdim])/2+1.78  )^.5
	
	return(list(
	"Theta.last"=Theta.last,
	"U.next"=U.next,
	"V.next"=V.next,
	"lambda.lasso"=lambda.lasso,
	"lambda.shrink"=lambda.shrink,
	"D.trunc"=D.trunc,
	"D.post"=D.post,
	"Theta.mode"=Theta.mode,
	"svd0"=svd0, 
	"Dtau"=Dtau,
	"D.ols"=svd0$d
	))
}



	expit<-function(x) exp(x)/(1+exp(x))
	
	################################
################################
##Three functions in here:
#### 1) test.gamma for finding the ML estimate for the grand lambda
#### 2) make.Z for converting a matrix to the Z-scale
#### 3) update.UDV for updating the ideal point estimates
#### *) Small additonal ones at the end 
################################
################################

################################
################################
##1) Finding gamma
#test.gamma.pois_gibbs<-function(gamma.try,Theta.last.0=Theta.last[row.type=="count",],votes.mat.0=votes.mat[row.type=="count",],emp.cdf.0,cutoff.seq=NULL){
test.gamma.pois_gibbs<-function(gamma.try,Theta.last.0,votes.mat.0,emp.cdf.0,cutoff.seq=NULL){
	gamma.try[gamma.try>2]<- 2
	gamma.try[gamma.try< -2]<- -2
	gamma.try<-exp(gamma.try)
	votes.mat<-votes.mat.0
	taus.try<-NULL
	count.seq<-cutoff.seq
	if(length(cutoff.seq)==0) count.seq<-seq(-1,max(votes.mat)+2,1)
	taus.try<-count.seq*0

	analytic.cdf<-count.seq
	a.int<-qnorm(mean(votes.mat==0))
	taus.try[count.seq>=0]<-(a.int+gamma.try[1]*count.seq[count.seq>=0]^(gamma.try[2]))
	taus.try[count.seq<0]<- -Inf
	taus.try<-sort(taus.try)
	taus.try[1]<--Inf	
	find.pnorm<-function(x){
		a<-x[1]
		b<-x[2]
		coefs<-c( -1.82517672, 0.51283415, -0.81377290, -0.02699400, -0.49642787, -0.33379312, -0.24176661, 0.03776971)
		x<-c(1, a, b, a^2,b^2, log(abs(a-b)),log(abs(a-b))^2,a*b)
		(sum(x*coefs))
	}
	lik<-((pnorm(taus.try[votes.mat+2 ]-Theta.last.0)-pnorm(taus.try[votes.mat+1 ]-Theta.last.0)) )
	log.lik<-log(lik)
	
	which.zero<-which(lik==0)
	a0<-taus.try[votes.mat[which.zero]+2]-Theta.last.0[which.zero]
	b0<-taus.try[votes.mat[which.zero]+1]-Theta.last.0[which.zero]
	log.lik[which.zero]<-apply(cbind(a0,b0),1,find.pnorm)
	log.lik[is.infinite(lik)&lik>0]<-max(log.lik[is.finite(lik)])
	log.lik[is.infinite(lik)&lik<0]<-min(log.lik[is.finite(lik)])

	thresh<-min(-1e30, min(log.lik[is.finite(log.lik)],na.rm=TRUE))
	log.lik[log.lik<thresh]<-thresh
	data.dev<--2*sum(log.lik,na.rm=TRUE)
	dev.out<-data.dev+sum(log(gamma.try)^2)
	
	
	return(list("deviance"=dev.out,"tau"=taus.try,"data.dev"=data.dev,"probvec"=log.lik))

}	



################################
################################
##2) Converting a matrix to the Z scale
	make.Z_gibbs<-function(
			Theta.last.0=Theta.last, 
			votes.mat.0=votes.mat,row.type.0=row.type, 
			n0=n, k0=k, params=NULL, iter.curr=0,empir=NULL,cutoff.seq.0=NULL,missing.mat.0=NULL,lambda.lasso,proposal.sd,scale.sd, max.optim,step.size,maxdim.0,tau.ord.0,test.gamma.pois.0
			){
		
		tau.ord<-tau.ord.0
		Theta.last<-Theta.last.0;votes.mat<-votes.mat.0;
		row.type<-row.type.0; n<-n0; k<-k0
		cutoff.seq<-cutoff.seq.0
		sigma<-1
		row.type<-row.type.0; votes.mat<-votes.mat.0
		Z.next<-matrix(NA,nrow=n,ncol=k)
		missing.mat<-missing.mat.0
		i.gibbs<-iter.curr
		maxdim<-maxdim.0
	test.gamma.pois<-test.gamma.pois.0

		#print(i.gibbs)
		#print(iter.curr)
	
	##Initialize deviance and summary stats
	word.dev<-vote.dev<-0
	wordprob.mat<-voteprob.mat<-NULL
	
			if(sum(row.type=="bin")>2){
				
			Z.next[row.type=="bin",][votes.mat[row.type=="bin",]==1]<-rtruncnorm(sum(votes.mat[row.type=="bin",]==1), mean=sigma^.5*Theta.last[row.type=="bin",][votes.mat[row.type=="bin",]==1], a=0, sd=sigma^.5)
		Z.next[row.type=="bin",][votes.mat[row.type=="bin",]==0]<-rtruncnorm(sum(votes.mat[row.type=="bin",]==0), mean=sigma^.5*Theta.last[row.type=="bin",][votes.mat[row.type=="bin",]==0], b=0 , sd=sigma^.5)
		Z.next[row.type=="bin",][missing.mat[row.type=="bin",]==1]<-rnorm(sum(missing.mat[row.type=="bin",]==1), mean=sigma^.5*Theta.last[row.type=="bin",][missing.mat[row.type=="bin",]==1], sd=sigma^.5)
		
		
		pars.max<-1
		accept.out<-prob.accept<-1
		Z.one<-Z.next[row.type=="bin",][votes.mat[row.type=="bin",]==1]
		Z.zero<-Z.next[row.type=="bin",][votes.mat[row.type=="bin",]==0]
		Z.one[abs(Z.one)>38]<-38*sign(Z.one[abs(Z.one)>38])
		Z.zero[abs(Z.zero)>38]<-38*sign(Z.zero[abs(Z.zero)>38])
		
		
		prob.vec<-as.vector(c(pnorm(-Z.zero,log=TRUE ),pnorm(Z.one,log=TRUE )))
		voteprob.mat<-data.frame("bin",exp(prob.vec),prob.vec,prob.vec^2)
		
		names(voteprob.mat)<-c("type","p","lp","psq")
		
		vote.dev<--2*sum(prob.vec,na.rm=T)
		
			
				}

	

		if(sum(row.type=="count")>2){
		#begin type = "pois"
		num.sample.count<-sum(row.type=="count")*k
		accept.out<-prob.accept<-NA

		dev.est<-function(x) test.gamma.pois_gibbs(x,votes.mat.0=votes.mat[row.type=="count",], Theta.last.0=Theta.last[row.type=="count",],cutoff.seq=cutoff.seq)$de
		dev.est1<-function(x) test.gamma.pois_gibbs(c(x,params[2]),votes.mat.0=votes.mat[row.type=="count",], Theta.last.0=Theta.last[row.type=="count",],cutoff.seq=cutoff.seq)$de
		dev.est2<-function(x) test.gamma.pois_gibbs(c(params[1],x),votes.mat.0=votes.mat[row.type=="count",], Theta.last.0=Theta.last[row.type=="count",],cutoff.seq=cutoff.seq)$de
	
	
	#range.opt<-c(params-1,params+1)
	grad.est<-function(x,del=.001){
		grad1<-(dev.est1(x[1]+del)-dev.est1(x[1]-del))/(4*del)
		grad2<-(dev.est2(x[2]+del)-dev.est2(x[2]-del))/(4*del)
		grad.all<-c(grad1,grad2)
		grad.all[!is.finite(grad.all)]<-0
		grad.all
		}
	grad.est.1<-function(x,del=.001){
		#grad1<-(dev.est1(x[1]+del)-dev.est1(x[1]-del))/(4*del)
		#grad2<-0#(dev.est2(x[2]+del)-dev.est2(x[2]-del))/(4*del)
		#c(grad1,grad2)
		grad1<-dev.est1(x[1]+del)-2*dev.est(x)+dev.est1(x[1]-del)
		grad1<-grad1/(2*del)
		grad1[!is.finite(grad1)]<-0
		grad1
		}
	grad.est.2<-function(x,del=.001){
		#grad1<-0#(dev.est1(x[1]+del)-dev.est1(x[1]-del))/(4*del)
		#grad2<-(dev.est2(x[2]+del)-dev.est2(x[2]-del))/(4*del)
		#grad2[!is.finite(grad2)]<-0
		grad2<-dev.est2(x[2]+del)-2*dev.est(x)+dev.est2(x[2]-del)
		grad2<-grad2/(2*del)
		#c(grad1,grad2)
		grad2
		}
	grad.est.12<-function(x,del=.001){
		#grad1<-0#(dev.est1(x[1]+del)-dev.est1(x[1]-del))/(4*del)
		#grad12<-(dev.est(x+del)-dev.est(x-del))/(8*del)
		grad12<-dev.est(x+del)-dev.est(x+del*c(1,-1))-dev.est(x+del*c(-1,1))+dev.est(x-del)
		grad12<-grad12/(4*del)
		grad12[!is.finite(grad12)]<-0
		grad12
		#c(grad1,grad2)
		}


			##Hamiltonian step
		##Code adapted from Neal, Handbook of MCMC
		gamma.opt<-"Hamiltonian Step"
		step.size<-min(step.size,2)
		
		#print("step size")
		#print(step.size)
		
		q0<-q<-params
		p0<-p<-rnorm(length(q),0,1)
		current_U<-dev.est(q0)/2

		num.HMC<-10
		for(i.HMC in 1:(num.HMC)){
		if(i.HMC%%3==1){
		hess.mat<-matrix(NA,nrow=2,ncol=2)
		hess.mat[1,1]<-grad.est.1(q)
		hess.mat[1,2]<-hess.mat[2,1]<-grad.est.12(q)
		hess.mat[2,2]<-grad.est.2(q)
		hess.mat[!is.finite(hess.mat)]<-0
		hess.mat<-.5*(hess.mat+t(hess.mat))
		epsilon<-ginv(hess.mat)*step.size/10
		}
		q<-q+epsilon%*%p
		p<-p-epsilon%*%grad.est(q)
		if(abs(dev.est(q)/2 - current_U)>10000) {revert.q<-TRUE;q<-q0;break}
		revert.q<-FALSE
		}
		q<-q+epsilon%*%p
		p<-p-epsilon%*%grad.est(q)/2
		p<- -p
		
		if(revert.q) q<-q0
		
		current_K<-sum(p0^2)/2
		proposed_U<-dev.est(q)/2
		proposed_K<-sum(p^2)/2
		proposal.dev<-proposed_U+proposed_K
		current.dev<-current_U+current_K
		gamma.cand<-q
		#print("Log Probs")
		#print(c(proposal.dev,current.dev))
		if(exp(current.dev-proposal.dev)>runif(1)){
			##Do if accept
			pars.max<-gamma.next<-(gamma.cand)
			#print("accept")
			accept.out<-prob.accept<-1
			} else {
			##Do if reject
			pars.max<-gamma.next<-(params)
			#print("reject")
			accept.out<-prob.accept<-0
			}
		
		#print("Proposal sd")
		#print(proposal.sd)
		#taus<-test.gamma.pois_gibbs(gamma.next,votes.mat.0=votes.mat[row.type=="count",], Theta.last.0=Theta.last[row.type=="count",],cutoff.seq=cutoff.seq)$tau
		pois.obj<-test.gamma.pois(gamma.next,votes.mat.0=votes.mat[row.type=="count",], Theta.last.0=Theta.last[row.type=="count",],cutoff.seq=cutoff.seq)
		taus<-pois.obj$tau
		#print(pois.obj$data)
		word.dev<-pois.obj$data
		prob.vec<-as.vector(pois.obj$probvec)
		
		wordprob.mat<-data.frame("count",exp(prob.vec),prob.vec,prob.vec^2)
		names(wordprob.mat)<-c("type","p","lp","psq")
		


		taus<-sort(taus)
		taus[1]<--Inf
		#print("Range of tau")
		#print(range.opt)
		#print(gamma.opt)
		#print(range(taus))
		#print("Beta parameters")
		#print(pars.max)
		
		Z.next.temp<-rtruncnorm(num.sample.count, mean=Theta.last[row.type=="count",], a=taus[votes.mat[row.type=="count",]+1 ],b=taus[votes.mat[row.type=="count",]+2 ] )

		
		Z.next[row.type=="count",]<-matrix(Z.next.temp,nrow=sum(row.type=="count"))



		}
	dev.all<-word.dev+vote.dev


	waicmat.all<-rbind(wordprob.mat,voteprob.mat)

	return(list("Z.next"=Z.next,"params"=(pars.max),"accept"=accept.out,"prob"=prob.accept,"proposal.sd"=proposal.sd,"step.size"=step.size,"tau.ord"=taus,"dev"=c(word.dev,vote.dev),"waicmat"=waicmat.all))
	}##Closes out make.Z function


################################
################################
##2) Converting a matrix to the Z scale
	update_UDV_gibbs<-function(
	Z.next.0=Z.next,
	k0=k, n0=n, lambda.lasso.0=lambda.lasso,lambda.shrink.0=lambda.shrink,
	Dtau.0=Dtau,
	votes.mat.0=votes.mat, iter.curr=0,row.type.0,missing.mat.0=missing.mat,maxdim.0,
	V.last,ratio.bin.0, prior.0
	){

	missing.mat<-missing.mat.0
	Dtau<-Dtau.0;
	Z.next<-Z.next.0;
	votes.mat<-votes.mat.0;
	n<-n0; k<-k0; 
	lambda.lasso<-lambda.lasso.0
	row.type <- row.type.0
	maxdim<-maxdim.0
	prior<-prior.0

	#Declare some vectors
	sigma<-1
	ones.r<-rep(1,k0)
	ones.c<-rep(1,n0)

	#Update intercepts
	mu.r<-rowMeans(Z.next)#-ones.c%*%t(mu.c)-Theta.last.0+mu.grand)
	mu.r<-mu.r*n/(n+1)+rnorm(length(mu.r),sd=1/k)

	mu.c<-colMeans(Z.next)#-mu.r%*%t(ones.r)-Theta.last.0+mu.grand)
	mu.c<-mu.c*k/(k+1)+rnorm(length(mu.c),sd=1/n)
	mu.grand<-mean(Z.next)
	mu.grand<-mu.grand*(n*k)/(n*k+1)+rnorm(1,sd=1/(n*k))
	
	mean.mat<-ones.c%*%t(mu.c)+mu.r%*%t(ones.r)-mu.grand
	
	if(length(unique(row.type))>1){
		is.bin<-sum(row.type=="bin")>0
		is.count<-sum(row.type=="count")>0
		is.ord<-sum(row.type=="ord")>0
		
		#print("Two Means Being Used")
		mean.c.mat<-matrix(NA,nrow=n,ncol=k)
		if(is.bin) mu.c.bin<-colMeans(Z.next[row.type=="bin",])
		if(is.count) mu.c.count<-colMeans(Z.next[row.type=="count",])
		if(is.ord) mu.c.ord<-colMeans(Z.next[row.type=="ord",])

		if(is.bin) mu.c.bin<-mu.c.bin*k/(k+1)+rnorm(length(mu.c.bin),sd=1/sum(row.type=="bin"))
		if(is.count) mu.c.count<-mu.c.count*k/(k+1)+rnorm(length(mu.c.count),sd=1/sum(row.type=="count"))
		if(is.ord) mu.c.ord<-mu.c.ord*k/(k+1)+rnorm(length(mu.c.ord),sd=1/sum(row.type=="ord"))

		if(is.bin) mean.c.mat[row.type=="bin",]<-ones.c[row.type=="bin"]%*%t(mu.c.bin)
		if(is.count) mean.c.mat[row.type=="count",]<-ones.c[row.type=="count"]%*%t(mu.c.count)
		if(is.ord) mean.c.mat[row.type=="ord",]<-ones.c[row.type=="ord"]%*%t(mu.c.ord)

	
		mean.grand.mat<-matrix(NA,nrow=n,ncol=k)
		if(is.bin) mean.grand.mat[row.type=="bin",]<-mean(Z.next[row.type=="bin",])+rnorm(1,sd=1/(sum(row.type=="bin")*k))
		if(is.count) mean.grand.mat[row.type=="count",]<-mean(Z.next[row.type=="count",])+rnorm(1,sd=1/(sum(row.type=="count")*k))
		if(is.ord) mean.grand.mat[row.type=="ord",]<-mean(Z.next[row.type=="ord",])* (sum(row.type=="count")*k)/(sum(row.type=="ord")*k+1)+rnorm(1,sd=1/(sum(row.type=="ord")*k))


	mean.mat<-mean.c.mat+mu.r%*%t(ones.r)-mean.grand.mat
	
	}

	Z.starstar<-svd.mat<- Z.next- mean.mat



	

	#Take svd, give each column an sd of 1 (rather than norm of 1)
	num.zeroes<-1#colMeans(votes.mat!=0)
	svd.mat.0<-svd.mat
	svd.mat.0<-apply(svd.mat.0,2,FUN=function(x) x-mean(x))
	drop.cols<-colMeans(svd.mat.0^2)^.5<1e-2
	drop.cols[!is.finite(drop.cols)]<-TRUE
	if(sum(drop.cols)>0) svd.mat[,drop.cols]<-rnorm(nrow(svd.mat)*sum(drop.cols))/3
	
	save(svd.mat,file="svd.mat")
	
	svd.mat[is.na(svd.mat)]<-0
	svd.mat.dum<-svd.mat
	
	wts.dum<-rep(1,nrow(svd.mat))
	wts.dum[row.type=="bin"]<-1/sum(1-missing.mat[row.type=="bin",])^.5*ratio.bin.0
	wts.dum[row.type=="count"]<-1/sum(1-missing.mat[row.type=="count",])^.5
	wts.dum[row.type=="ord"]<-1/sum(1-missing.mat[row.type=="ord",])^.5


	wts.dum<-wts.dum/mean(wts.dum)

	svd.dum<-svd(svd.mat*wts.dum,nu=maxdim,nv=maxdim)
	#svd.dum<-rsvd(svd.mat*wts.dum,k=maxdim)

	#svd.dum<-irlba(svd.mat*wts.dum,nu=maxdim,nv=maxdim,V=V.last)
	svd.dum$u[!is.finite(svd.dum$u)]<-0
	svd.dum$d[!is.finite(svd.dum$d)]<-0
	svd.dum$v[!is.finite(svd.dum$v)]<-0


	svd.dum$d<-(t(svd.dum$u)%*%svd.mat%*%svd.dum$v)
	which.rows<-which(rowMeans(svd.dum$d^2)^.5<1e-4)
	which.cols<-which(colMeans(svd.dum$d^2)^.5<1e-4)
	svd.dum$d[which.rows,]<-rnorm(length(svd.dum$d[which.rows,]),sd=.001)
	svd.dum$d[which.cols,]<-rnorm(length(svd.dum$d[which.cols,]),sd=.001)
	svd2<-svd(svd.dum$d,nu=maxdim,nv=maxdim)
	#svd2<-rsvd(svd.dum$d,nu=maxdim,nv=maxdim,k=maxdim)

	#print(dim(svd.dum$d))
	#svd2<-irlba(svd.dum$d,nu=maxdim,nv=maxdim)
	svd2$u[is.na(svd2$u)|is.infinite(svd2$u)]<-0
	svd2$d[is.na(svd2$d)|is.infinite(svd2$d)]<-0
	svd2$v[is.na(svd2$v)|is.infinite(svd2$v)]<-0
	
	svd.dum$u<-svd.dum$u%*%(svd2$u)
	svd.dum$v<-svd.dum$v%*%(svd2$v)

	svd.dum$d<-svd2$d
	
	svd0<-svd.dum
	
	svd0$v<-t(t(svd0$u)%*%svd.mat.0)
	svd0$v<-apply(svd0$v,2,FUN=function(x) x/sum(x^2)^.5)
	svd0$d<-diag(t(svd0$u)%*%svd.mat.0%*%svd0$v)
	sort.ord<-sort(svd0$d,ind=T,decreasing=T)$ix
	svd0$u<-svd0$u[,sort.ord]
	svd0$v<-svd0$v[,sort.ord]
	svd0$d<-svd0$d[sort.ord]
	svd0$u<-svd0$u*(n-1)^.5
	svd0$v<-svd0$v*(k-1)^.5
	svd0$d<-svd0$d*((n-1)*(k-1))^-.5
	Theta.last.0<-svd.mat
	Theta.last<-Theta.last.0+(ones.c%*%t(mu.c)+mu.r%*%t(ones.r)-mu.grand)

	#Update d; follows from Blasso and DvD
	Y.tilde<-as.vector(svd.mat)
	if(n>length(Dtau)) Dtau[(length(Dtau)+1):n]<-1
	A<- (n*k)*diag(n)+diag(as.vector(Dtau^(-1)))
	#if(n>k) for(i.A in (k+1):n) A[i.A,i.A]<-(Dtau^-1)[i.A]*0+1e20
	#gA<-ginv(A)
	gA<-A*0+NA
	gA[1:maxdim,1:maxdim]<-ginv(A[1:maxdim,1:maxdim])
	gA[is.na(gA)]<-0
	XprimeY<-sapply(1:maxdim, FUN=function(i, svd2=svd0, Z.use=svd.mat) sum(Z.use*(svd2$u[,i]%*%t(svd2$v[,i]))))
	if(length(XprimeY)<dim(gA)[1]) XprimeY[(length(XprimeY)+1) :dim(gA)[1]]<-0
	D.post.mean<- as.vector(gA%*%XprimeY)
	D.post.var.2<-gA

##Sample D and reconstruct theta, putting intercepts back in
	D.post<-rep(NA,length(D.post.mean))
	D.post[1:maxdim]<-as.vector(mvrnorm(1, mu=D.post.mean[1:maxdim], D.post.var.2[1:maxdim,1:maxdim] ) )/sigma^.5
	D.post[is.na(D.post)]<-0
	abs.D.post<-abs(D.post)

##Calculate MAP and mean estimate
	D.trunc<-pmax(svd0$d-lambda.lasso.0,0)
	
	U.last<-svd0$u
	V.last<-svd0$v

##Update U, V
	prior.var<-.25	
	U.last<-t(t(U.last)*(D.post[1:maxdim]^2/(D.post[1:maxdim]^2+1/(4*k))))
	ran.mat<-sapply(1/(D.post[1:maxdim]^2*(n-1)+.25),FUN=function(x) rnorm(nrow(U.last),mean=0,sd=x^.5))
	U.last<-U.last+ran.mat
	V.last<-t(t(V.last)*D.post[1:maxdim]^2/(D.post[1:maxdim]^2+1/(4*n)))
	ran.mat<-sapply(1/(D.post[1:maxdim]^2*(k-1)+.25 ),FUN=function(x) rnorm(nrow(V.last),mean=0,sd=x^.5))
	V.last<-V.last+ran.mat

	U.next<-U.last
	V.next<-V.last
	
	#Uncommenting the next two lines eliminates the U and V
	#Gibbs update.
	#U.next<-svd0$u
	#V.next<-svd0$v

	Theta.last.0<-U.next%*%diag(D.post[1:maxdim])%*%t(V.next)
	Theta.last<-Theta.last.0+mean.mat
	Theta.last[row.type=="bin",]<-Theta.last[row.type=="bin"]-mean(Theta.last[row.type=="bin",][missing.mat[row.type=="bin"]==0])+ qnorm(mean(votes.mat[row.type=="bin",][missing.mat[row.type=="bin"]==0]))
	#Not sure what the next three lines did.
	#chisq.ran<-rchisq(1,n*k)#+rchisq(maxdim,1)
	#sigma<-(sum((Z.next-Theta.last)^2))/sum(chisq.ran)
	#Theta.last<-Theta.last
	Theta.mode<- U.next%*%diag(D.trunc[1:maxdim])%*%t(V.next)+mean.mat
	##Update muprime, invTau2, lambda.lasso
	muprime<-(abs(lambda.lasso*sqrt(sigma)/(abs.D.post)))
	invTau2<-sapply(1:maxdim, FUN=function(i) rinv.gaussian(1, muprime[i], (lambda.lasso^2 ) ) )
	Dtau<-abs(1/invTau2)
	lambda.lasso<-rgamma(1, shape=maxdim+prior[1]  , rate=sum(Dtau[1:maxdim])/2+prior[2]  )^.5
	lambda.shrink<-((maxdim+prior[1]-1)/(sum(Dtau[1:maxdim])/2+prior[2]))^.5
	
	
	return(list(
	"Theta.last"=Theta.last,
	"U.next"=U.next,
	"V.next"=V.next,
	"lambda.lasso"=lambda.lasso,
	"lambda.shrink"=lambda.shrink,
	"D.trunc"=D.trunc,
	"D.post"=D.post,
	"Theta.mode"=Theta.mode,
	"svd0"=svd0, 
	"Dtau"=Dtau,
	"D.ols"=svd0$d
	))
}



	expit<-function(x) exp(x)/(1+exp(x))
	
	
###################################################################
###################################################################
## Function for generating the criterion
###################################################################
###################################################################
crit.function<-function(x,longrun=TRUE){
	null.use<-FALSE
	ratio.use<-10^c(3,1,.75,.5,.25,0,-.25,-.5,-.75,-1,-3,10)
	Sen.use<-105:112
	pars.all<-cbind(rep(ratio.use,times=length(Sen.use)),rep(Sen.use,each=length(ratio.use)))

	CongressNum<-pars.all[x,2]
	ratio.use<-pars.all[x,1]
	if(ratio.use==10^10) null.use<-TRUE

	file.out<-paste(CongressNum,ratio.use,sep="_")
	votes.name<-paste("votemat_",CongressNum,sep="")
	load(votes.name)	
	names.file<-paste("dwnom",CongressNum,".csv",sep="")
	names.mat<-read.csv(names.file)
	dtm.file<-paste("dtm_",CongressNum,sep="")
	load(dtm.file)


	set.seed(1)

	#votes.mat.0<-votes.mat
	#votes.mat.0[votes.mat.0==0]<-.5
	#votes.mat.0[votes.mat.0==6]<-0
	#votes.mat.0[votes.mat.0>1]<-0.5
	#votes.mat<-votes.mat.0

	dtm3<-rbind(t(votes.mat),t(dtm2))
	missing.mat.use<-rbind(t(votes.mat==0.5),t(dtm2)*0)*1

	scale.text<-NULL


	set.seed(1)
	if(longrun){
	scale.text<-criterion(dtm3,missing.mat.use, maxdim=20,max.optim=1000,gibbs=2000,burnin=2000, ratio=ratio.use)
	}else{
	scale.text<-criterion(dtm3,missing.mat.use,max.optim=10,gibbs=10,burnin=20, ratio=ratio.use)
	}
 disc.stat<-function(obj){
	d.temp<-rowMeans(obj$dim.mean)
 	tempstat<-0
 	for(i in 1:length(d.temp)) {
	      u.curr<-obj$coldims.all[,i]
		tempstat<-tempstat+sum(outer(u.curr,u.curr,"-")^2)*d.temp[i]
			}
				tempstat
 }
 
wmat.all<-scale.text$waic.out
wmat.words<-wmat.all[wmat.all[,1]=="count",]
wmat.vote<-wmat.all[wmat.all[,1]=="bin",]

waic.vote<- -2 * sum(log(wmat.vote[, 2]), na.rm = TRUE) + 
            2 * (sum(wmat.vote[, 4] - wmat.vote[, 3]^2))

waic.words<- -2 * sum(log(wmat.words[, 2]), na.rm = TRUE) + 
            2 * (sum(wmat.words[, 4] - wmat.words[, 3]^2))

stat.out<-c(disc.stat(scale.text),scale.text$WAIC,waic.vote,waic.words)

#save(stat.out,file=file.out)
return(c(stat.out,scale.text$WAIC))

}



	
